function [ results_fnames ] = osl_run_first_level_continuous_state_vector_seed( oat )

% [ results_fnames ] = osl_run_first_level_continuous( oat )
%
% takes in an OAT, which needs to be setup by calling oat=osl_setup_oat(S), struct
% and runs first level continuous time GLM or seed based connectivity
%
% This function should normally be called using osl_run_oat(oat);
%
% MWW 2012

global OSLDIR;

use_class_probs=0;

first_level=oat.first_level;
source_recon_name='source_recon';

% check if user has asked not to do GLM
if(~isfield(first_level,'doGLM')),
    first_level.doGLM = 1;
end

if(first_level.time_average)
    error('Time averaging over whole window is not appropriate on continuous data');   
end;

if ~isempty(first_level.tf_downsample_factor)
    first_level.tf_downsample_factor=1;
    warning('Turning off downsampling - not implemented');
end;

% set the downsample factor.  If doing a tf, read tf_downsample_factor.  If
% doing time domain, read time_downsample_factor
switch first_level.tf_method
    case 'none' % use time_downsample_factor
        error('Continuous data time-wise GLM only handles time-freq analysis');
        
        if ~isempty(first_level.time_downsample_factor)
            ds_factor = 1./first_level.time_downsample_factor;
        else
            ds_factor = [];
        end
    otherwise  % use tf_downsample_factor
        if ~isempty(first_level.tf_downsample_factor)
            ds_factor = 1./first_level.tf_downsample_factor;
        else
            ds_factor = [];
        end
end % switch first_level.tf_method

if(~isfield(first_level,'connectivity_seed_mni_coord')),
    first_level.connectivity_seed_mni_coord=[];
end;

Ntrials   = 1;
   
clear stats;

for subi_todo=1:length(first_level.sessions_to_do),   
        
    sub=first_level.sessions_to_do(subi_todo);
    
    %sub/length(first_level.submatfiles),
    
    %% load in the beamformer result for this session
    source_recon_results_fname = oat.source_recon.results_fnames{sub};
    source_recon_results=osl_load_oat_results(oat,source_recon_results_fname);
    
    % get D and update path in case OAT dir has been moved
    D = osl_get_oat_sensordata( source_recon_results )
    
    %% setup things that are common to all sessions
    if(subi_todo==1),
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %% set up mask indices from mask/coords
        if strcmp(source_recon_results.recon_method,'none'), % sensor space analysis
            
            disp('Working in sensor space');
  
            %%%%%%%%%%%%%%%
            %% setup channels
            modality='MEG';
            chanindmeg = strmatch(modality, D.chantype); % need to even include bad channels in output so that it is the same for all sessions
       
            first_level_results.mask_indices_in_source_recon=1:length(chanindmeg);
                   

        else,
    
            if(isfield(first_level,'mni_coords'))    
                                
                current_level_mni_coord = first_level.mni_coords;
                
                % load in lower level coords
                lower_level_mni_coord=source_recon_results.mni_coord;               
                %lower_level_mask_fname  = [oat.source_recon.dirname '/' source_recon_name '_mask'];
                %[ lower_level_mni_coord, ~ ] = osl_mnimask2mnicoords(lower_level_mask_fname);

                [~, iA, mask_indices]=intersect(current_level_mni_coord, lower_level_mni_coord,'rows');

                % sort indices to be in the order of the lower level mni_coords:
                [~, gg]=sort(iA);
                indices_in_lower_level=mask_indices(gg);

                first_level_results.mask_indices_in_source_recon = indices_in_lower_level;
                first_level_results.mni_coord = lower_level_mni_coord(indices_in_lower_level,:);
                
                % the mni coordinates from the specified list
                str=['Desired MNI coordinates ' ];
                for vox=1:size(current_level_mni_coord,1),
                    str=[str ', [' num2str(current_level_mni_coord(vox,:)) ']'];
                end;
                disp(str);
                
                % the mni coordinates from the lower level mask
                str=['Nearest dipoles from lower level at MNI coordinates ' ];
                for vox=1:length(indices_in_lower_level),
                    str=[str ', [' num2str(lower_level_mni_coord(indices_in_lower_level(vox),:)) ']'];
                end;
                disp(str);
            
            else
                % setup std space brain
                first_level_results.gridstep=source_recon_results.gridstep;

                S=[];
                S.lower_level_mask_fname=[oat.source_recon.dirname '/' source_recon_name '_mask'];
                S.current_level_mask_fname=[oat.source_recon.dirname '/' oat.first_level.name '_mask'];
                S.current_level=first_level;
                S.lower_level_mni_coord=source_recon_results.mni_coord;
                S.lower_level_gridstep=source_recon_results.gridstep;

                [first_level_results.mask_indices_in_source_recon, first_level_results.mni_coord]=setup_mask_indices(S);
                lower_level_mni_coord=source_recon_results.mni_coord
                clear S;
            end;
        end;       
    
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %% setup HMM 
        if strcmp(source_recon_results.recon_method,'none'),
            NK=1;
        else
            NK=numel(source_recon_results.BF.inverse.W.MEG); 
        end;
    
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %% set up potential time frequency bases in tf_settings, note that any time or
        % frequency averaging will be done AFTER the time-frequency decomposition
        
        % set the frequency range
        freq_range = first_level.tf_freq_range; % range used for stat testing
        if isempty(freq_range)
            % if not set as a first_level option, default to the range used
            % for source_recon
            freq_range=source_recon_results.source_recon.freq_range;
            
            if isempty(freq_range)
                error('No frequency range has been specified. Set oat.source_recon.freq_range or oat.first_level.tf_freq_range');
            end % isempty(freq_range)
        end % if isempty(freq_range)
        
        source_recon_time_indices=find(source_recon_results.samples2use);
        source_recon_times = D.time(source_recon_time_indices);
        
        if isempty(first_level.time_range)
            first_level.time_range=source_recon_results.woi;
        end
        
        tf_settings = [];
        tf_settings.tf_method               = first_level.tf_method;
        tf_settings.tf_logtransform         = first_level.tf_logtransform;
        tf_settings.tf_hilbert_freq_res     = first_level.tf_hilbert_freq_res;
        tf_settings.tf_freq_range           = freq_range; % set above
        tf_settings.tf_num_freqs            = first_level.tf_num_freqs;
        tf_settings.raw_times               = source_recon_times;
        tf_settings.ds_factor               = 1; % do not downsample inside transform - instead we will do this later after the moving average low pass filtering
        tf_settings.tf_morlet_factor        = first_level.tf_morlet_factor;
        tf_settings.tf_multitaper_taper                = first_level.tf_multitaper_taper;
        tf_settings.tf_multitaper_twin                 = first_level.tf_multitaper_twin;
        tf_settings.tf_multitaper_freqsmooth           = first_level.tf_multitaper_freqsmooth;
        tf_settings.tf_calc_amplitude=0;
        
        out = osl_tf_transform( tf_settings ); % passed without a data matrix, osl_tf_transform just gets container sizes
        
        if isfield(out, 'tf_morlet_basis')
            first_level.tf_morlet_basis = out.tf_morlet_basis;
        end
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %% setup time and frequency windows/indices for using later
          
        % setup result time dimension        
        first_level_time_indices  = ( out.tf_times >= first_level.time_range(1) ) & ( out.tf_times <= first_level.time_range(2) );
        first_level_results.times=1;       
        if first_level.time_average
            disp('Note time averaging irrelevant and so is ignored for continuous data time-wise GLM');
        end
        
        % setup result freq dimension
        first_level_results.frequencies = out.tf_freqs;
        if(first_level.freq_average),
            first_level_results.frequencies=mean(first_level_results.frequencies);
        end;
       
        % setup result spatial dimension
        Nvoxels_out=length(first_level_results.mask_indices_in_source_recon);
         
        clear out

        %%%%%%%%%%%%%%%%
        %% setup seed for seed-based connectivity analysis       
        if(isempty(first_level.connectivity_seed_mni_coord)),
           error('need to specify seed');
        end;
        
        if(first_level.tf_num_freqs>1)
            error('Continuous data seed-based connectivity only handles single freq bands at a time');
        end;

        mni_coords=lower_level_mni_coord;
        
        disp(['Using for seed: MNI coordinate ' num2str(first_level.connectivity_seed_mni_coord)]);

        % find the nearest index of the beamformed voxels to the specified mni_coord
        dists=(sqrt(sum((mni_coords-repmat(first_level.connectivity_seed_mni_coord,size(mni_coords,1),1)).^2,2)));

        [dist,first_level_results.seed_index]=min(dists);

        disp(['Nearest dipole for seed at MNI coordinate ' num2str(mni_coords(first_level_results.seed_index,:))]);                                              

        first_level.contrast=[];
        first_level.contrast{1}=[1]';

        %%%%%%%%%%%%%%%%
        %% results containers
        ntpts=1;
        nfreqs=length(first_level_results.frequencies);

        first_level_results.pseudo_zstat_var=zeros([Nvoxels_out,1]);
        
    end; % end of if(sub==1)
    
    disp(['Reconstruct time courses and computing stats for dataset ' source_recon_results.fname]);    

    %%%%%%%%%%%%%%%
    %% setup channels
    modality='MEG';
    chanindmeg = strmatch(modality, D.chantype);
    chanind = setdiff(chanindmeg, D.badchannels);
    if isempty(chanind)
        error(['No good ' modality ' channels were found.']);
    end
    if strcmp(source_recon_results.recon_method,'none'), % sensor space analysis
        % for sensor space - also include bad sensors sensors 
        first_level_results.mask_indices_in_source_recon=1:length(chanindmeg);
        chanind = chanindmeg;
    end;    
    
    %%%%%%%%%%%%%%%
    %% load in sensor data
    first_level_results.D_sensor_data=D;
    first_level_results.chanind=chanind;
    
    triallist=1;
    
    sensor_data=D(chanind, source_recon_time_indices, 1);   
    
    classchanind=find(strcmp(D.chanlabels,'Class'));
    if(isempty(classchanind)),
        error(['No ''Class'' chanlabel in: ' D.fname]);
    else
        for si=1:NK,        
            class_samples_inds{si} = (D(classchanind,source_recon_time_indices, 1)==si);  

            state_tinds=(class_samples_inds{si}==1);                                                

            total_time=(source_recon_times(end)-source_recon_times(1))*size(state_tinds,3);
            state_time=total_time*sum(sum(squeeze(state_tinds)))/prod(size(squeeze(state_tinds)));
            txt=['State ' num2str(si) ' is active for ' num2str(state_time) 'secs'];

            if(state_time<10),
                warning([txt]);
            else
                disp(txt);
            end;
        
        end; 
    end;

    %%%%%%%%%%%%%%%%%%%%%%
    %% recon the seed for doing orthogonalisation  
    if(isempty(first_level.connectivity_seed_mni_coord)),
        error('No seed specified');
    end;
        
    % recon seed time course
    S2=[];
    S2.source_recon_results=source_recon_results;
    S2.chanind=chanind;
    S2.class_samples_inds=class_samples_inds;
    S2.voxind=first_level_results.seed_index;
    S2.use_class_probs=0;
    S2.sensor_data=sensor_data;
    S2.reduce_rank=2;

    [ dat_seed wnorms_seed wnorm_seed wnorms_nai_seed wnorm_nai_seed, weights_seed ] = get_voxel_recon_timecourse_vector( S2 );

    %dat_seed=dat_seed';
    %wnorms_nai_seed=wnorms_nai_seed';

    dat_seed=dat_seed./sqrt(mean(squash(wnorms_nai_seed))); % trials x tpts x dipole_rank

    conn.dat_seed=permute(dat_seed(:,first_level_time_indices,:),[2 3 1]);
    
    clear S2;
    
    %%%%%%%%%%%%%%%
    %% Orthogonalise sensor space data wrt seed:
    sensor_data_zerolagcorr=zeros(size(sensor_data));
    if(first_level.connectivity_seed_regress_zerolag),

        do_seed_regress_zerolag_statewise=0;
        
        % orthoganalise the target voxel WRT the seed voxel (doing
        % it this way round means that the TF transform and
        % windowing does not need to be redone for the seed voxel)
        if(NK>1 && do_seed_regress_zerolag_statewise),  
            
            for si=1:NK,
                tinds=(class_samples_inds{si}==1);
                
                for tri=1:size(sensor_data,3),
                
                    dat2=permute(dat_seed(tri,tinds,:),[2,3,1]);

                    % iterate over sensors
                    for indind=1:size(sensor_data,1),
                        dat1=permute(sensor_data(indind,tinds,tri),[2,1,3]);

                        if(indind~=first_level_results.seed_index)

                            % iterate over dipole directions
                            for dd=1:size(dat_seed,3),
                                dat2=permute(dat_seed(tri,tinds,dd),[2,3,1]);
                                dat1=dat1-dat2*(dat2\dat1);
                            end;

                        end;

                        sensor_data_zerolagcorr(indind,tinds,tri)=dat1;
                    end;
                
                end;
            end;
        else,
                     
            for tri=1:size(sensor_data,3),
                
                dat2=permute(dat_seed(tri,:,:),[2,3,1]);

                % iterate over sensors
                for indind=1:size(sensor_data,1),
                    dat1=permute(sensor_data(indind,:,tri),[2,1,3]);
                                       
                    if(indind~=first_level_results.seed_index)
                        
                        % iterate over dipole directions
                        for dd=1:size(dat_seed,3),
                            dat2=permute(dat_seed(tri,:,dd),[2,3,1]);
                            dat1=dat1-dat2*(dat2\dat1);
                        end;
                        
                    end;
                    
                    sensor_data_zerolagcorr(indind,:,tri)=dat1;
                end;
                
            end;
        end;        

    end;        
    
    %%%%%%%%%%%%%%%
    %% Do T-F transform in sensor space
    %% Transform the data - either downsampling in time domain, or a TF transform
    
    disp('Doing T-F transform in sensor space...');

    tf_settings = [];
    tf_settings.tf_method               = first_level.tf_method;
    tf_settings.tf_logtransform         = first_level.tf_logtransform;
    tf_settings.tf_hilbert_freq_res     = first_level.tf_hilbert_freq_res;
    tf_settings.tf_freq_range           = freq_range; % set above
    tf_settings.tf_num_freqs            = first_level.tf_num_freqs;
    tf_settings.raw_times               = source_recon_times;
    tf_settings.ds_factor               = 1; % do not downsample inside transform - instead we will do this later after the moving average low pass filtering
    tf_settings.tf_morlet_factor        = first_level.tf_morlet_factor;
    tf_settings.tf_multitaper_taper                = first_level.tf_multitaper_taper;
    tf_settings.tf_multitaper_twin                 = first_level.tf_multitaper_twin;
    tf_settings.tf_multitaper_freqsmooth           = first_level.tf_multitaper_freqsmooth;
    tf_settings.tf_calc_amplitude=0;

    if isfield(first_level,'tf_morlet_basis')
        tf_settings.tf_morlet_basis = first_level.tf_morlet_basis;
    end

    for indind=1:size(sensor_data_zerolagcorr,1),

        for tri=1:size(sensor_data_zerolagcorr,3),
            out = osl_tf_transform( tf_settings , sensor_data_zerolagcorr(indind,:,tri) ); % do the transformation
            sensor_data_zerolagcorr_tf(indind,:,tri)=out.dattf;
        end;
    end;
        
    %%%%%%%%%%%%%%%%
    %% Now do windowing over time                    
    sensor_data_zerolagcorr=sensor_data_zerolagcorr(:,first_level_time_indices,:);
    sensor_data_zerolagcorr_tf=sensor_data_zerolagcorr_tf(:,first_level_time_indices,:);
    
    %%%%%%%%%%%%%%
    %% Get TF for seed data 
    S2=[];
    S2.source_recon_results=source_recon_results;
    S2.chanind=chanind;
    S2.class_samples_inds=class_samples_inds;
    S2.voxind=first_level_results.seed_index;
    S2.use_class_probs=0;
    S2.sensor_data=sensor_data_zerolagcorr_tf;
    S2.reduce_rank=2;

    [ dat_seed_tf wnorms_seed wnorm_seed wnorms_nai_seed wnorm_nai_seed, weights_seed ] = get_voxel_recon_timecourse_vector( S2 );
        
    conn.dat_seed_tf=abs(permute(dat_seed_tf(1,:,:),[2 3 1]));
    clear S2;
    
    %%%%%%%%%%%%%%%%
    %% do everything from now on one voxel at a time (for the sake of
    %% RAM) in the first level mask space                
    disp('Looping over voxels...');

    ft_progress('init', 'etf');

    % debug:
    %[indfind, mni_coord]=nearest_vec(source_recon_results.mni_coord,[-32 -26 50]);
            
    if any(strcmp(oat.first_level.connectivity_fc_methods,'he')) || any(strcmp(oat.first_level.connectivity_fc_methods,'plv')) || any(strcmp(oat.first_level.connectivity_fc_methods,'synchro_hmm')) || any(strcmp(oat.first_level.connectivity_fc_methods,'hmm')),
    for indind=1:Nvoxels_out, % indexes brain space      
        
        ft_progress(indind/Nvoxels_out);                                   

        S2=[];
        S2.source_recon_results=source_recon_results;
        S2.chanind=chanind;
        S2.class_samples_inds=class_samples_inds;
        S2.voxind=first_level_results.mask_indices_in_source_recon(indind);
        S2.sensor_data=sensor_data_zerolagcorr;
        S2.use_class_probs=0;
        S2.reduce_rank=2;

        [ dat wnorms wnorm wnorms_nai wnorm_nai weights] = get_voxel_recon_timecourse_vector( S2 );

        dat=dat./sqrt(mean(squash(wnorms_nai)));
           
        S2.sensor_data=sensor_data_zerolagcorr_tf;
        [ dat_tf wnorms wnorm wnorms_nai wnorm_nai weights] = get_voxel_recon_timecourse_vector( S2 );

        dat_tf=dat_tf./sqrt(mean(squash(wnorms_nai))); 
        
        conn.dat_tf=permute(dat_tf(1,:,:),[2 3 1]);
        conn.dat=permute(dat(1,:,:),[2 3 1]);
    
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
        %% now do the stats

        clear state_tinds;

        %%%%%%
        % on vec:

        fsample=1/(source_recon_times(2)-source_recon_times(1));   
        tres=1/fsample;

        bandcentre=mean(first_level.tf_freq_range);

        if any(strcmp(oat.first_level.connectivity_fc_methods,'synchro_hmm')),
            %%
            bandcentre=mean(first_level.tf_freq_range);
            S=[];

            S.tres=tres;
            S.data_high=(conn.dat_seed);
            S.data_low=(conn.dat);
            S.NK=4;
            S.force_zero_means=1;
            S.num_hmm_starts=1;
            S.do_plots=0;
            S.hmm_pca_dim=10;                
            span=1/bandcentre;
            S.M=round(span/tres);
            S.deltasecs=span/S.M;

            %S.deltasecs=1/(S.M*bandcentre);

            res  = synchro_hmm( S );

            kl=zeros(NK);

            hmm=res.hmm;
            Apca=res.Apca;
            NK=S.NK;
            for ii=1:NK,     
                for jj=1:NK,
                    if(ii~=jj),
                        ci=hmm.state(ii).Cov;
                        cj=hmm.state(jj).Cov;
                        kl(ii,jj)=trace(pinv(cj)*ci)-log(det(ci)/det(cj))-size(cj,1);
                    end;
                end;
            end; 

            first_level_results.pw_oasm_embedded(indind,si)=sum(squash(kl))/(NK^2-NK);
            %%
        end;

        if any(strcmp(oat.first_level.connectivity_fc_methods,'hmm')),
            %%
            bandcentre=mean(first_level.tf_freq_range);
            S=[];

            S.tres=tres;
            S.data_high=(conn.dat_seed);
            S.data_low=(conn.dat);
            S.NK=4;
            S.force_zero_means=1;
            S.num_hmm_starts=1;
            S.do_plots=0;
            S.hmm_pca_dim=4;                
            S.deltasecs=tres;
            S.M=1;
            %S.deltasecs=1/(S.M*bandcentre);

            res  = synchro_hmm( S );

            kl=zeros(NK);

            hmm=res.hmm;
            Apca=res.Apca;
            NK=S.NK;
            for ii=1:NK,     
                for jj=1:NK,
                    if(ii~=jj),
                        ci=hmm.state(ii).Cov;
                        cj=hmm.state(jj).Cov;
                        kl(ii,jj)=trace(pinv(cj)*ci)-log(det(ci)/det(cj))-size(cj,1);
                    end;
                end;
            end; 

            first_level_results.pw_oasm(indind,si)=sum(squash(kl))/(NK^2-NK);

            %%
        end;

        %%%%%%
        %% search for dir
        w=weights{si}';   

        if(size(w,2)>1)
        else
            oat.first_level.connectivity_seed_search_direction=0;
        end;

        tmp=w'*w;   

        [u, ~] = svd(tmp,'econ');                              
        yeta1 = (u(:,1));

        if(~oat.first_level.connectivity_seed_search_direction),
            yeta2 = zeros(size(yeta1));
        else
            yeta2 = (u(:,2));
        end;

        w=weights_seed{si}';                
        tmp=w'*w;   

        [u, ~] = svd(tmp,'econ');                                  
        xeta1 = (u(:,1));

        if(~oat.first_level.connectivity_seed_search_direction),
            xeta2 = zeros(size(xeta1));
        else
            xeta2 = (u(:,2));
        end;

        if(oat.first_level.connectivity_seed_search_direction),
            as=0:0.4:pi;                
            bs=0:0.4:pi;               
        else
            as=0;
            bs=0;
        end;
        
        if(indind>0),

            plv=zeros(length(as),length(bs));
            he=zeros(length(as),length(bs));

            clear cope varcope;
            for a=1:length(as),

                % now do windowing over time                    
                x=conn.dat_seed_tf;
                xin=(x*(xeta1*cos(as(a))+xeta2*sin(as(a))));

                % now compute HE:
                x=(sqrt(xin.*conj(xin)));

                % moving average 
                if ~isempty(first_level.time_moving_av_win_size)

                    windowSize = first_level.time_moving_av_win_size*fsample;

                    x=moving(x,round(windowSize));

                end

                x=normalise(x);
                pinvxtx=pinv(x'*x);
                pinvx=pinv(x);

                x_phase=angle(xin);

                xraw=conn.dat_seed;
                xraw=(xraw*(xeta1*cos(as(a))+xeta2*sin(as(a))));

                for b=1:length(bs),

                    if any(strcmp(oat.first_level.connectivity_fc_methods,'he')) || any(strcmp(oat.first_level.connectivity_fc_methods,'plv')),
                        
                        % now do windowing over time                    
                        y=conn.dat_tf;
                        yin=(y*(yeta1*cos(bs(b))+yeta2*sin(bs(b))));

                        % now compute HE:
                        y=(sqrt(yin.*conj(yin)));

                        % moving average 
                        if ~isempty(first_level.time_moving_av_win_size)

                            fsample=1/(source_recon_times(2)-source_recon_times(1));        
                            windowSize = first_level.time_moving_av_win_size*fsample;

                            y=moving(y,round(windowSize));

                        end

                        y=normalise(y);

                        y_phase=angle(yin);

                        [cope, varcope]=glm_fast_for_meg(y,x,pinvxtx,pinvx,[1]);

                        he(a,b)=cope/sqrt(varcope);

                        plv(a,b) = fisher_xform(abs(mean(exp(sqrt(-1)*(y_phase-x_phase))))); % Measuring Phase Synchrony in Brain Signals. Jean-Philippe Lachaux et al, HBM 1999
                    end;

                    if any(strcmp(oat.first_level.connectivity_fc_methods,'synchro')),
                        M=10;
                        deltasecs=1/(M*bandcentre);
                        delta=ceil(deltasecs/tres);

                        yraw=conn.dat;
                        yraw=(yraw*(yeta1*cos(bs(b))+yeta2*sin(bs(b))));
                        gs(a,b)=synchro([normalise(yraw) normalise(xraw)],M,delta,10,10);
                    end;
                end;

            end;

            tmp=he;

            [i,j]=find(tmp==max(max(tmp)));                
            first_level_results.he_angles_seed(indind,si)=as(i);                
            first_level_results.hes(indind,si)=max(max(tmp));

            tmp=plv;

            [i,j]=find(tmp==max(max(tmp)));                
            first_level_results.plv_angles_seed(indind,si)=as(i);                
            first_level_results.plvs(indind,si)=max(max(tmp));

            if any(strcmp(oat.first_level.connectivity_fc_methods,'synchro')),

                tmp=gs;

                [i,j]=find(tmp==max(max(tmp)));                
                first_level_results.gs_angles_seed(indind,si)=as(i);                
                first_level_results.gss(indind,si)=max(max(tmp));
            
            end;
            
        end;           
                  
    end; % for indind = 1:Nvoxels_out
    end; % if any(strcmp(oat.first_level.connectivity_fc_methods
    
    if any(strcmp(oat.first_level.connectivity_fc_methods,'hmm_kl')),
    
        disp('Doing HMM...');           
                 
        NK=6;
        hmm_num_starts=1;
        hmm_kl_use_source_seed=0;
        hmm_pca_dim=40;
        use_seed=1;
        
        disp(['Inferring ' num2str(NK) ' HMM states']);            
        
        Sh=[];                

        if(hmm_kl_use_source_seed)
            sensor_data_zerolagcorr_tmp=permute(sensor_data_zerolagcorr,[2 1 3]);
            %sensor_data_zerolagcorr_tmp=permute(sensor_data,[2 1 3]);
            sd=[size(sensor_data_zerolagcorr_tmp) 1 1]; 
            sensor_data_zerolagcorr_tmp=reshape(sensor_data_zerolagcorr_tmp,[sd(1), sd(2)*sd(3)]);
            Sh.data=sensor_data_zerolagcorr_tmp;
            Sh.data=[Sh.data, conn.dat_seed];
            
            if(hmm_pca_dim>0)
                disp('Doing pre-HMM PCA...');
                [allsvd,Apca]=pca(Sh.data,hmm_pca_dim);
                pinvApca=pinv(Apca);              
            else
                error('hmm_pca_dim unspecified');
            end;
            
            Sh.data=(pinvApca*Sh.data')';    
        else,
            sensor_data_tmp=permute(sensor_data,[2 1 3]);
            %sensor_data_zerolagcorr_tmp=permute(sensor_data,[2 1 3]);
            sd=[size(sensor_data_tmp) 1 1]; 
            sensor_data_tmp=reshape(sensor_data_tmp,[sd(1), sd(2)*sd(3)]);
            Sh.data=sensor_data_tmp;
            
            if(hmm_pca_dim>0)
                disp('Doing pre-HMM PCA...');
                [allsvd,Apca]=pca(Sh.data,hmm_pca_dim);
                pinvApca=pinv(Apca);              
            else
                error('hmm_pca_dim unspecified');
            end;
            
            Sh.data=(pinvApca*Sh.data')';    
        end;
        
        % add seed onto the end:
        
        Sh.data=normalise(Sh.data);
        
        Sh.NK=NK;

        Sh.num_starts=hmm_num_starts;

        [ hmm, block ]=run_multistart_hmm(Sh);                
        
        NK=hmm.K;       
        
        Shh=[];
        Shh.NK=NK;
        Shh.hmm=hmm;        
        Shh.block=block;
        Shh.tres=(source_recon_times(2)-source_recon_times(1));   
        Shh.Apca=Apca;
        
        plot_hmm( Shh );
  
        cv=Apca*hmm.state(1).Cov*Apca';
        
        if(hmm_kl_use_source_seed)
            weights_seed=zeros(2,size(cv,1));
            weights_seed(:,size(cv,1)-1:end)=eye(2);    
        else
            weights_seed=source_recon_results.BF.inverse.W.MEG{1}{first_level_results.seed_index}; % needs to be changed to get working with HMM beamformer
            tmp=weights_seed*weights_seed';   
            [u, ~] = svd(tmp,'econ');                              
            weights_seed = u(:,1:2)'*weights_seed;
            weights_seed=weights_seed/sqrt(trace(weights_seed*weights_seed'));
            weights_seed=[weights_seed];
        end;
                
        if(numel(source_recon_results.BF.inverse.W.MEG)>1),
            error('Not implemented for HMM beamformer');
        end;
        
        if(use_seed)
            kl=zeros(Nvoxels_out,NK,NK);
        
            ft_progress('init', 'etf');
            for voxind=1:Nvoxels_out;

                ft_progress(voxind/Nvoxels_out);                                   

                clear cv_pair;
                for kk=1:NK,

                    cv=Apca*hmm.state(kk).Cov*Apca';

                    weights=source_recon_results.BF.inverse.W.MEG{1}{voxind}; % needs to be changed to get working with HMM beamformer
                    tmp=weights*weights';   
                    [u, ~] = svd(tmp,'econ');                              
                    weights = u(:,1:2)'*weights;
                    weights=weights/sqrt(trace(weights*weights'));

                    if(hmm_kl_use_source_seed)
                        weights=[weights, zeros(2,2)];
                    end;

                    weights_pair=[weights', weights_seed']';

                    cv_pair{kk}=weights_pair*cv*weights_pair';   

                end;            

               for jj=1:NK,
                    cj=cv_pair{jj};
                    pcj=inv(cj);
                    ldcj_mNj=logdet(cj,'chol')-size(cj,1);

                    for ii=jj+1:NK,   
                    %for jj=1:NK,
                    %    if(ii~=jj),    
                        ci=cv_pair{ii};

                        %  kl(voxind,ii,jj)=trace(pinv(cj)*ci)-log(det(ci)/det(cj))-size(cj,1);
                        kl(voxind,ii,jj)=trace(pcj*ci)-logdet(ci,'chol')+ldcj_mNj;  
                    %end;
                    end;
                end;   
                                                                        
            end;  
            
            ft_progress('close');
        else,
                    
            % precompute weights
            for voxind=1:Nvoxels_out;
                
                weights=source_recon_results.BF.inverse.W.MEG{1}{voxind}; % needs to be changed to get working with HMM beamformer
                tmp=weights*weights';   
                [u, ~] = svd(tmp,'econ');                              
                weights = u(:,1:2)'*weights;
                weights=weights/sqrt(trace(weights*weights'));
                weights_stored{voxind}=weights;
            end;
            
            % precompute state cvs:
            for kk=1:NK,
                cvs_stored{kk}=Apca*hmm.state(kk).Cov*Apca';
            end;
            
            kl=zeros(Nvoxels_out,Nvoxels_out,NK,NK);
        
            ft_progress('init', 'etf');
            
            cv_pair=zeros(size(weights_stored{1},1)*2,size(weights_stored{1},1)*2,NK);

            for voxind_seed=1:Nvoxels_out;

                ft_progress(voxind_seed/Nvoxels_out);  

                for voxind=1:Nvoxels_out;

                    if(voxind_seed~=voxind),                                                                     

                        clear cv_pair;
                        for kk=1:NK,

                            weights_pair=[weights_stored{voxind}', weights_stored{voxind_seed}']';
                            cv_pair(:,:,kk)=weights_pair*cvs_stored{kk}*weights_pair';   

                        end;            


                        for jj=1:NK,
                            cj=cv_pair(:,:,jj);
                            pcj=inv(cj);
                            ldcj_mNj=logdet(cj,'chol')-size(cj,1);

                            for ii=jj+1:NK,     
                                ci=cv_pair(:,:,ii);

                                %  kl(voxind,ii,jj)=trace(pinv(cj)*ci)-log(det(ci)/det(cj))-size(cj,1);
                                kl(voxind_seed,voxind,ii,jj)=trace(pcj*ci)-logdet(ci,'chol')+ldcj_mNj;  
                                
                            end;
                        end;        
                    end;
                end;
            end;
           
            ft_progress('close');
        end;
                
        %% node OASM
        % cd ~; load kl;
        
        first_level_results.kl_edge_oasm=sum(sum(kl,3),4);
        
        first_level_results.kl_node_oasm=sum(sum(sum(kl,2),3),4);
   
        first_level_results.klm_edge_oasm=max(max(kl,[],3),[],4);
        
        first_level_results.klm_node_oasm=max(max(max(kl,[],2),[],3),[],4);
      
        first_level_results.kl_edge_index=(1:size(first_level_results.kl_edge_oasm,1))';
    end;
    
    if (first_level.doGLM && ~isempty(first_level.connectivity_seed_mni_coord)),
        first_level_results.cope_by_state(first_level_results.seed_index,:,:,:)=1;
        first_level_results.stdcope_by_state(first_level_results.seed_index,:,:,:)=1e16;
        first_level_results.cope(first_level_results.seed_index,:,:,:)=1;
        first_level_results.stdcope(first_level_results.seed_index,:,:,:)=1e16;
    end;
                
    first_level_results.source_recon_results_fname=source_recon_results_fname;
    first_level_results.first_level=first_level;
    first_level_results.source_recon=source_recon_results.source_recon;
    first_level_results.level=1;
    %first_level_results.source_recon_results_fname=source_recon_results.fname;
    first_level_results.name=first_level.name;
    first_level_results.recon_method=source_recon_results.recon_method;
    first_level_results.session_name=source_recon_results.session_name;
          
    first_level_results.fname=[ first_level_results.session_name '_' first_level_results.name ];
    disp(['Saving statistics in file ' oat.source_recon.dirname '/' first_level_results.fname]);
    osl_save_oat_results( oat, first_level_results );
    results_fnames{sub}=first_level_results.fname;

    if (strcmp(source_recon_results.recon_method,'none')),
    else
        disp(['To create niftii files from this use a call to osl_save_nii_stats']);
    end;
    
end;

function y=fisher_xform(x)

y=0.5*log((1+x)/(1-x));
    

